\newcommand{\Oalpha}{\Omega}
\newcommand{\Oalphae}{\Omega^e}
\newcommand{\Oalphaobs}{\Omega^*}
\newcommand{\Galpha}{\Gamma}
\newcommand{\bhh}{\mathit{h}}
\newcommand{\bgg}{\mathit{g}}
\newcommand{\Galphagg}{\Gamma_{\bgg}}
\newcommand{\Galphahh}{\Gamma_{\bhh}}
\newcommand{\Galphag}{\Gamma_{\bg}}
\newcommand{\Galphah}{\Gamma_{\bh}}
\newcommand{\Galphatgd}{\Gamma_{g,d}}
\newcommand{\Galphathd}{\Gamma_{h,d}}
\newcommand{\Galphath}{\Gamma_{\bh}}

\input{gbmacros.tex}

\chapter{Biophysical Model of the Heart}
\label{chap:model}

In this chapter we describe the anatomical structure of the human heart and describe the aspects that are taken into consideration while constructing the
biophysical model. The model does not attempt to be a generic all encompassing model of the heart, but attempts to factor key aspects of cardiac dynamic that are
important for the successful estimation of myocardial motion. The use of this model to constrain the problem of cardiac motion estimation shall be explained in
detail in Chapter \ref{chap:inverse}. The heart is a complicated system, with an electro-mechanical system responsible for the activation and contraction of the
heart muscles. Modeling of cardiac anatomy, electrophysiology and mechanics is an active research field and a comprehensive review of the field can be found in
\cite {sachse04}.

Muscle fiber orientations need to be considered while modeling cardiac electromechanics. The diffusion properties of the muscle fibers play an important role in
the propagation of cardiac activation current. Similarly the force generated by the muscles is primarily along the fiber direction. Therefore knowledge of the
diffusion properties within the myocytes or at least the myocyte orientations is very important for modeling purposes, especially if the patient specific models
are desired. Diffusion Tensor MR Imaging (DTI)  \cite{scollan98, tseng99} is a promising new imaging modality that allows the diffusion tensors within tissues to
be quantified. Additionally, the principal direction of the diffusion tensor aligns with myocyte orientations \cite{sachse01, hsu01}. However, it is not possible
to obtain cardiac DTI {\em in vivo} currently, and as a result most modeling approaches use synthetic fiber orientations \cite{sachse04}. Synthetic fiber
orientations are in general very simple and capture only the basic trends in the orientation. Instead of using synthetic models for myocyte orientations we
estimate patient specific myocyte orientation using {ex vivo} DTI obtained from a human volunteer. The process of warping the {\em ex vivo} DT to patient
specific DT shall be described in detail in Section \ref{sec:warpDTI}.

Cardiac myocytes on being excited contract along the myocyte orientation. This contraction causes the atria and ventricles to contract, in turn causing blood to
ciculate within the body. We model the myocardium, blood and surrounding coelomic cavity as linear elastic materials and the active forces within the myocardium
caused as a result of the myocyte contraction. The linear elastodynamic model is described in Section \ref{sec:linElas} and the formulation of the active forces
is described in Section \ref{sec:forces}. Finally we present experimental results validating the various models described in this Chapter in Section
\ref{sec:modelResults}. We start by giving a brief overview of cardiac anatomy.

\section{Anatomical Structure of the Heart}
\label{sec:anatomy}

Modeling of cardiac anatomy, electrophysiology, and mechanics is very important for understanding the complicated interactions that take place between different
anatomical structures. This knowledge helps us to understand the mechanisms of heart failure, and can help devise ways to prevent and cure such pathologies. A
large number of cardiac pathologies occur because of problems with the electro-mechanical system within the heart. Consequently, a lot of active research is
being carried out in this field.  A comprehensive review of the field can be found in \cite{sachse04}.

The walls of the heart are composed of cardiac muscle, called myocardium. It consists of four compartments: the right and left atria and ventricles. The heart is
oriented so that the anterior aspect is the right ventricle while the posterior aspect shows the left atrium. The left ventricular free wall and the septum are
much thicker than the right ventricular wall. This is logical since the left ventricle pumps blood to the systemic circulation, where the pressure is
considerably higher than for the pulmonary circulation, which arises from right ventricular outflow. Since a muscle fiber can contract only in one direction, the
heart structure is complex, to succeed at pumping the blood. Anatomically, to achieve this, the muscle walls of the ventricles and the atria are composed of a
single helically folded muscular structure as can be seen in Figure \ref{grays_1}. The cardiac muscle fibers are divided into four groups \cite{gray18}: Two
groups of fibers wind around the outside of both ventricles. Beneath these fibers a third group winds around both ventricles. Beneath these fibers a fourth group
winds only around the left ventricle.  

As a result of this complicated geometry, muscle fiber orientations need to be considered while modeling cardiac electro-mechanics. At present the prevailing hypothesis is that the diffusion properties of the muscle fibers play an important role in the propagation speed of cardiac activation current as the diffusion tensor is proportional to the local conductivity, and also in the estimation of the orientation of forces generated by the muscles which are along the fiber direction. Therefore, knowledge of the diffusion tensor or at least the fiber orientations is very important for modeling purposes, especially if patient specific models are desired. However, {\em in-vivo} diffusion tensor imaging of the heart is not a practically realizable application at present, and as a result most modeling approaches use synthetic data for the fiber orientations. A common approach is to vary the elevation angle between the fiber and the short axis plane between $+90^{\circ}$ and $-90^{\circ}$ from the endocardium to the epicardium \cite{sachse04, sermesant04}. These models of fiber orientations capture only the overall trends in orientation, and thus are not sufficient for patient specific modeling.

%Since the myocyte orientations are responsible for both the electrical conduction and the mechanical contraction of the myocardium, it is important to consider the fiber orientations while modeling the heart. 

\begin{figure}[!hbtp]
\centering
\includegraphics[width=0.4\textwidth]{images/grays_3}
\includegraphics[width=0.34\textwidth]{images/grays_2} 
\caption{Fiber orientations in the Human Heart showing the helical structure of the muscles (from \cite{gray18})}
\label{grays_1}
\end{figure} 

\input{fiber}

\section{Mechanical Modeling}

We shall now describe the linear elastic formulation of the mechanical model of the heart. We model the heart as a linear elastic solid occupying a bounded region $\omega$ and assume that the normalized fiber orientations, $\bn$, at the resting phase of the heart (end of diastole) are available.

We can write the forward problem as:
%%
\begin{equation}
\label{eq:model} M\ddot{u}(t) + C\dot{u}(t) + Ku(t) + F(t) = 0\quad
 t \in (0,1).
\end{equation} 

%%

Using a Ritz-Galerkin formulation, with basis
 functions $\phi$, the expressions for $M$ and $K$
 are given by $M_{ij}=\int {\mathrm I} (\phi_i \phi_j)$,
 $K=\int (\lambda + \mu) \Grad\phi_i \otimes \Grad \phi_j + \mu
 {\mathrm I} (\Grad\phi_i \cdot \Grad\phi_j)$, and $C=\alpha M + \beta K$, with $\alpha$
 and $\beta$ viscous-damping parameters.  Here, $\otimes$ is the outer
 vector product, $\lambda$ and $\mu$ are the Lam\'{e} constants, and
 ${\mathrm I}$ is the 3D identity matrix. Equation \eqref{eq:model} is
 derived by the Navier linear elastodynamics equation
 \cite{gurtin-81}. 

To solve (\ref{eq:model}) we embed $\omega$ in a regular domain $\Omega$. We use trilinear finite elements to discretize (\ref{eq:model}), and piecewise constant functions for $\lambda,\nu$.  The Poisson ratio $\nu({\bf x})$ varies between $0$ for a fully compressible material to $0.5$ for a fully incompressible material. When considering soft tissue deformations, the value of the Poisson ratio given for many tissue classes borders on the limit of incompressibility\footnote{As $\nu$ approaches 0.5, commonly used displacement-based finite element implementations suffer from the so-called locking effect.  We use underintegration for the $\nabla \cdot {\bf u}$ term in (\ref{eq:model}); see \cite{Hughes87} for details.}. Gladilin~\cite{Gladilin2003} studies the sensitivity of the Poisson ratio on the displacement field obtained via an incompressible formulation and a compressible formulation. The results indicate that Poisson ratio does not affect the solution significantly.

\begin{figure}
  \centering  
  \includegraphics[width=0.5\textwidth]{images/mat_prop}
  \caption{The material properties and boundary conditions for the cardiac model. We set the stress, $\sigma = 0$, at the boundary $\Gamma$. Different Lam\'{e}  parameters are selected for the myocardium ($\lambda_1, \mu_1$), blood ($\lambda_2, \mu_2)$, the lungs ($\lambda_3, \mu_3)$ and for bone ($\lambda_4, \mu_4)$ }
 \label{matprop}
\end{figure}	

We consider the input data as defined by the MR/DTI image to define the regular domain $\Omega$. The boundary conditions are exactly prescribed on the six faces of the cuboid, $\Gamma$, are simply set to have zero stress, i.e., $\sigma = 0$. One important issue is the accurate integration of the elements that overlap in the inhomogeneity transition. For simplicity we have used voxelized material	properties even in the case of exact geometry information. Within the regular domain $\Omega$, different material properties are assigned based on the tissue type. We currently consider the myocardium ($\lambda_1, \mu_1$), blood ($\lambda_2, \mu_2$), the lungs ($\lambda_3, \mu_3$) and bone ($\lambda_4, \mu_4$). These are shown in Figure \ref{matprop}. The material properties for the different tissue types are listed in Table \ref{tab:matprop}.

\begin{table}[tbp]
	\begin{center}
		\begin{tabular}{|l|c|c|c|c|}
		\hline
		Tissue Type & myocardium & blood & lungs & bone \\
		\hline
		Young's Modulus (kPa) & 10 & 1  & 1 & 100 \\
		Poisson's ratio & 0.45   & 0.35 & 0.35 &  0.45 \\
		\hline
		\end{tabular}
	\end{center}
	\caption{The material properties used for the different tissue types used in our cardiac model.}
	\label{tab:matprop}
\end{table}


\input{fem}
\input{forces}

\section{Validation}
\label{sec:modelResults}

In this section we present results from the experiments that were conducted to validate the fiber estimation and the numerical correctness of the mechanical model.

\subsection {Validation of the Fiber Estimation Algorithm}

We used canine DTI datasets obtained from Center for Cardiovascular Bioinformatics and Modeling, Johns Hopkins University and acquired at the National Institue of Health, to validate the effectiveness of our diffusion tensor remapping algorithm. A total of $19$ canine subjects were scanned, of which $12$ subjects were normal, and $7$ had cardiac failure. The scans were performed {\em in vitro} after the hearts were harvested. Each heart was placed in an acrylic container filled with Fomblin, a perfluoropolyether (Ausimon, Thorofare, NJ). Fomblin has a low dielectric effect and minimal MR signal, thereby increasing contrast and eliminating unwanted susceptibility artifacts near the boundaries of the heart.  The long axis of the hearts was aligned with the $z$ axis of the scanner.  Images were acquired with a 4-element knee phased array coil on a 1.5 T GE CV/I MRI Scanner (GE, Medical System, Wausheka, WI) using an enhanced gradient system with 40 mT/m maximum gradient amplitude and a 150 T/m/s slew rate.

Since only the long axis of the heart was aligned with the $z$ axis of the scanner, we first need to correct for rotation about the $z$ axis. We picked one of the normal canine hearts as the template, and performed affine registration to warp the remaining subjects onto the template space. The diffusion tensors for these subjects were also rotated by the rotational component of the affine transform. We then perform the elastic registration using HAMMER to estimate the transformation that maps the template to the subject space. This transformation is then used to warp and reorient the diffusion tensors of template onto the subject. The quality of the mapping is measured by computing the angle between the principal direction of the mapped tensor and the principal direction of the ground truth obtained from the subject's DTI. This is shown in Figure \ref{fig:dot}.

The error in the fiber orientations is further demonstrated on one slice in Figure \ref{fig:pd_both}, by comparing the original fibers (in red) and the mapped fibers (in blue). Our method was able to successfully map the fibers for healthy as well as for failing canine hearts. The percentage of voxels where the error in the principal directions is less than $10^\circ$ is shown in Figure \ref{fig:graph}. We evaluate an error of less than $10^\circ$ since it is close to the average error obtained from DTI imaging \cite{scollan98} and histological measurements\cite{streeter73}. We also compare this with the error by using a synthetic model to estimate fiber orientations. The synthetic fiber orientations were produced by varying the elevation angle between the short axis plane and the fiber between $+90^\circ$ and $-90^\circ$ from the endocardium to the epicardium. Similar synthetic models have been used in \cite{sachse04, sermesant04}.

\begin{figure}
\begin{center}
\framebox{
%\includegraphics[width=0.4\textwidth]{images/dot_prod} 
\includegraphics[width=0.4\textwidth]{images/dot}}
\caption{The angle between the mapped principal direction with the actual principal direction, in degrees. The image is overlaid on a segmentation of the heart based on the fractional anisotropy } 
\label{fig:dot}
\end{center}  
\end{figure} 

The method was able to map the fibers accurately for both healthy as well as for failing hearts that had left bundle branch blocks. Although, the percentage of voxels having an error of less than $10^\circ$ was lower in the failing hearts as compared to the healthy ones, we observe that most of these errors are in the vicinity of the block as expected, and the fiber orientations away from the bundle branch block were similar to that observed in healthy patients. Therefore, our method should perform better for other pathologies like ventricular hypertropy and infarction, since it has been shown that the muscle fiber orientations are not affected significantly by hypertrophy and myocardial infarction\cite{walker05}. 

\begin{figure}
\begin{center}
%\includegraphics[width=0.4\textwidth]{images/dot_prod} 
\includegraphics[width=0.5\textwidth]{images/f1} 
\includegraphics[width=0.395\textwidth]{images/fbr3} 
\caption{The principal directions of the original DT are shown in blue. The mapped principal directions are shown in red. The glyphs are overlaid on a segmentation of the heart based on the fractional anisotropy of the mapped DT}
\label{fig:pd_both}
\end{center}  
\end{figure} 

\begin{figure}
\begin{center}
%\includegraphics[width=0.4\textwidth]{images/dot_prod} 
\includegraphics[width=0.8\textwidth]{images/graph} 
\caption{The percentage of voxels having less than $10^\circ$ error in the principal directions after mapping }
\label{fig:graph}
\end{center}  
\end{figure} 

\subsection{Synthetic Models for the validation of the forward solver}

In order to validate the numerical accuracy and correctness of the forward, we used two synthetic models for validation. The first is a gaussian distribution of the forces with homogeneous material properties. The material properties are uniform over the entire domain and we choose a material withPoisson's ration $\nu=0.45$ and Young's modulus $1$ kPa. There are no fiber orientations in this model, and we work with full forces (\ref{e:genericForce}). We define forces over this domain as,

\begin{equation}
f_x = f_y = f_z = sin(kt) e^{-(x^2+y^2+z^2)}.
\end{equation}


In order to assess the fiber-orientation parametrized model of the forces (\ref{e:fiberForce}), we use an ellipsoidal model of the left ventricle. The fiber orientations are generated by varying the elevation angle between the fiber and the short axis plane between $+60^{\circ}$ and $-60^{\circ}$ from the endocardium to the epicardium \cite{sachse04,sermesant04}. The model along with the fiber orientations is shown in Figure \ref{f:fiberModel}. For this model we selected a Poisson's ratio $\nu=0.45$ and a Young's modulus of 10 kPa for the myocardial tissue and 1 kPa for the surrounding tissue and ventricular cavity. Raleigh damping ($C =\alpha M + \beta K$) was used with parameters $\alpha=0$ and $beta=7.5\times10^{-4}$. In order to drive the forward model, we generated
forces by propagating a synthetic activation wave from the apex to the
base of the ventricles. Snapshots of this activation wave at different
phases of the cardiac cycle are shown in Figure \ref{f:activation}.

\begin{figure}
	\centering
		\subfloat[]{
		\includegraphics[width=0.4\textwidth]{figs/fiberModel2}}
		\subfloat[]{
		\includegraphics[width=0.4\textwidth]{figs/fiberModel1}}
	\caption{The ellipsoidal LV model shown with ({\bf a}) the fiber orientations, and ({\bf b}) the fibers obtained via fiber tracking.}
	\label{f:fiberModel}
\end{figure}

\begin {figure}[tbp]
\centering
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 0mm, clip,width=0.25\textwidth]{figs/act0}
	}
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 0mm, clip,width=0.25\textwidth]{figs/act1}
	}
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 0mm, clip,width=0.25\textwidth]{figs/act2}
	} \\
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 0mm, clip,width=0.25\textwidth]{figs/act3}
	}
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 0mm, clip, width=0.25\textwidth]{figs/act4}
	}
\caption{The activation function used in the
synthetic model of the heart. The activation function starts at the
apex and moves to the base of the ventricles.}
\label{f:activation}
\end{figure}



The octree is generated using a gaussian distribution of points for the first case and using the geometry for the synthetic LV model. The balanced octrees for both models are shown in Figure \ref{f:octrees}. In both cases the smallest element size was set to be the same size as the regular grid counterpart. The number of octants for different grid sizes for both the models is listed in Table \ref{tab:fwdRes}.

\begin{figure}
	\centering
		\subfloat[] {\includegraphics[width=0.4\textwidth]{figs/gaussianOct}}
		\subfloat[] {\includegraphics[width=0.4\textwidth]{figs/ellipsoidOct}}
	\caption{The octrees corresponding to ({\bf a}) the gaussian force model, and ({\bf b}) the ellipsoidal fiber force model.}
	\label{f:octrees}
\end{figure}

\subsection{Numerical Verification of the Forward Solver}

We use the gaussian model to study the convergence of the forward
solver by comparing the results obtained from the numerical solver
with the analytical solution. The solution was obtained for
discretizations corresponding to regular grids $32^3$, $64^3$, $128$, and
$256^3$. The number of time steps were equal to the number of elements (along a axis) in all cases.
The relative residual tolerance for the CG solver was set to $10^-16$ for the gaussian model. The time steps were again equal to the number of elements and the relative residual tolerance for the CG solver was set to $10^-8$ 
For both models the use of octrees for meshing reduced the number of elements significantly. We
observed $\mathcal{O}(h^2)$ convergence, where $h$ is the grid
size. These results, along with the number of processors used and the
time for the computation are tabulated in Table \ref{tab:fwdRes}.

\begin{table}
\scriptsize
\centering
		\begin{tabular}{|c|c|c|c|l|l|c|c|l|}\hline
		Problem &  Number of	&	\multicolumn{4}{c|}{Full Force} & \multicolumn{3}{c|}{Fiber Force} \\\cline{3-9}
			Size	&	Processors & \#Octants & \#Iters & Time &	Error &\#Octants & \#Iters & Time \\\hline
$32^3$ &	1  &	7386 & 5  & 87ms	& $7.63\times10^{-2}$	& 5958 & 19 & 291ms \\
$64^3$ &	8  &	52K & 6  & 126ms	& $1.91\times10^{-2}$	& 24K & 21 & 323ms  \\
$128^3$	 &	64  & 396K &	8 & 186ms	&	$4.77\times 10^{-3}$ & 100K & 26 & 383ms \\
$256^3$	 &	512 &	 2.8M & 11 & 244ms & $1.21\times 10^{-3}$ & 456K & 31 &	443ms \\
\hline	 
		\end{tabular}
\caption{Results from the validation of the forward solver for the two models. The regular grid size and the number of processors that the problem was run on are listed. We also list the number of octants in the resulting balanced Octree. The average number of KSP iterations and wall clock time per time step are listed. The relative error($\|\cdot\|_2$)are also provided. For the gaussian force, the relative error is computed with respect to the analytical solution.
\label{tab:fwdRes}}
\end{table}

 The second
model used was constructed from a MR image of a human heart. The
fiber orientation was obtained from ex-vivo Diffusion Tensor (DTI)
images of the heart. In order to drive the forward model, we generated
forces by propagating a synthetic activation wave from the apex to the
base of the ventricles. Snapshots of this activation wave at different
phases of the cardiac cycle are shown in Figure \ref{fig:activation}.
The fiber orientation, the myocardial forces and the resulting
deformations within the myocardium are shown at different slices and
time points in Figures \ref{fig:fibers} and \ref{fig:fibers1}.

\begin {figure}[tbp]
\centering
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 50mm, clip, width=0.3\textwidth]{images/fwd/act0}
	}
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 50mm, clip, width=0.3\textwidth]{images/fwd/act1}
	}
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 50mm, clip, width=0.3\textwidth]{images/fwd/act2}
	} \\
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 50mm, clip, width=0.3\textwidth]{images/fwd/act3}
	}
	\subfloat[]{
		\includegraphics[trim = 50mm 50mm 50mm 50mm, clip, width=0.3\textwidth]{images/fwd/act4}
	}
\caption{The activation function used in the
synthetic model of the heart. The activation function starts at the
apex and moves to the base of the ventricles.}
\label{fig:activation}
\end{figure}


\begin {figure}[tbp]
\centering
	\subfloat[]{
		\includegraphics[width=0.3\textwidth]{images/fwd/fibers}
	}
	\subfloat[]{
		\includegraphics[width=0.3\textwidth]{images/fwd/forces}
	}
	\subfloat[]{
		\includegraphics[width=0.3\textwidth]{images/fwd/deformation}
	}
\caption{(a) The orientation of the cardiac fibers,
(b) the forces developed within the myocardium as a result of the
contraction of the fibers, and (c) the resulting deformations within
the myocardium.}
\label{fig:fibers}
\end{figure}

\begin {figure}[tbp]
\centering \subfloat[]{
	\includegraphics[width=0.3\textwidth]{images/fwd/def_34} }
	\subfloat[]{
	\includegraphics[width=0.3\textwidth]{images/fwd/def_44} }
	\subfloat[]{
	\includegraphics[width=0.3\textwidth]{images/fwd/def_54} }
\caption{The deformations induced by the propagation of
the activation from the apex to the base.}
\label{fig:fibers1}
\end{figure}


\section{Conclusions}

In this chapter we introduced algorithms for estimating myocardial fiber orientations and described the linear elastic bio-mechanical model of the heart that we shall use for constraining the motion estimation problem. This is a test.
